function SPM = CIATstats(pathArray, ageArray)
    [pathstr, name, ext] = fileparts(pathArray{1});
    % Set P,I 
    P = pathArray;
    nScan = length(P);
    I = ones(nScan,4); %set the group appartenance
    for i=2:nScan
        I(i,1) = i-1;
        I(i,2) = 2;
    end
    VY    = spm_vol(char(P));
    H = [zeros(nScan,1) ones(nScan,1)];
    H(1,:) = [1 0];
    Hnames = [{'Group_1'} {'Group_2'}];
    B = ones(nScan,1);
    C = [];
    Cnames = {};
    G = ageArray - spm_meanby(ageArray,ones(nScan,1));
    g=zeros(nScan,1);
    for i = 1:nScan 
        g(i) = spm_global(VY(i)); 
    end
    gSF    = 50./g; %global mean set to 50
    xC(1)= struct('rc',ageArray,'rcname','Age','c',G,'cname','Age',...
                'iCC',1,'iCFI',1,'type',2,'cols',4, ...
                'descrip',{{'nuisance variable: Age';'used centered around overall mean'}});
    xC(2)= struct('rc',g,'rcname','global','c',{[]},'cname',{{}},...
            'iCC',0,'iCFI'	,0,'type',3,'cols',{[]},'descrip',...
            {{'global values: (used for proportional scaling)';...
            '("raw" unscaled globals shown)'}});
    for i = 1:nScan
        VY(i).pinfo(1:2,:) = VY(i).pinfo(1:2,:)*gSF(i);
    end
    %iGloNorm : see spm_spm_ui.m
    %Global scaling options
    % 8: grand mean scaling
    xGX = struct('iGXcalc',3,'sGXcalc',...
        'mean voxel value (within per image fullmean/8 mask)',...
        'rg',g,'iGMsca',8,	'sGMsca',...
        '(implicit in PropSca global normalisation)',	'GM',50,'gSF',gSF,...
        'iGC',12,'sGC','(redundant: not doing AnCova)','gc',{[]},...
        'iGloNorm',8,'sGloNorm','proportional scaling to 50 ');
    
    %-Construct full design matrix (X), parameter names and structure (xX)
    %===================================================================
    X      = [H C B G];
    tmp    = cumsum([size(H,2), size(C,2), size(B,2), size(G,2)]);
    xX     = struct(	'X',		X,...
        'iH',		[1:size(H,2)],...
        'iC',		[1:size(C,2)] + tmp(1),...
        'iB',		[1:size(B,2)] + tmp(2),...
        'iG',		[1:size(G,2)] + tmp(3),...
        'name',		{[Hnames Cnames {'\mu'} {'Age'}]},...
        'I',		I,...
        'sF',		{{'Subject' 'Group' '' ''}});
    
    % Non-sphericity correction (set xVi.var and .dep)
    %===================================================================
    xVi.I   = I;
    nL      = max(I);		                         % number of levels
    mL      = find(nL > 1);		                     % multilevel factors
    xVi.var = sparse(1,4);                           % unequal variances
    xVi.dep = sparse(1,4);                           % dependencies
    xVi     = spm_non_sphericity(xVi);
    
    mask_threshold = 0.8;
    xsM.Analysisthreshold = ['images thresholded at' num2str(mask_threshold) 'times global'];
    xsM.Explicit_masking = 'No';
    mask_threshold = 0.5;
    xM     = struct('T',0.+mask_threshold*sqrt(-1), 'TH',mask_threshold*50*ones(nScan,1),...
        'I',0, 'VM',{[]}, 'xs',xsM);
    tmp = [ {'2 condition, +0 covariate, +1 block, +1 nuisance'};...
            {'4 total, having 3 degrees of freedom'};...
            {'leaving 4 degrees of freedom from 7 images'} ];
    xsDes = struct(	'Design', 'Compare-populations: 1 scan/subject (AnCova)',...
        'Global_calculation','mean voxel value (within per image fullmean/8 mask)',...
        'Grand_mean_scaling','(implicit in PropSca global normalisation)',...
        'Global_normalisation',	'proportional scaling to 50',...
        'Parameters',{tmp});
            
    %-Assemble SPM structure
    %===================================================================
    SPM.xY.P	= P;			% filenames
    SPM.xY.VY	= VY;			% mapped data
    SPM.nscan	= nScan;        % scan number
    SPM.xX		= xX;			% design structure
    SPM.xC		= xC;			% covariate structure
    SPM.xGX		= xGX;			% global structure
    SPM.xVi		= xVi;			% non-sphericity structure
    SPM.xM		= xM;			% mask structure
    SPM.xsDes	= xsDes;		% description
    SPM.SPMid	= 'ICEM v1';	% version
    SPM.swd     = pathstr;      % SPM working directory
    
    %-Save SPM.mat
    %-------------------------------------------------------------------
    save(fullfile(pathstr,'SPM.mat'), 'SPM');